# robustness tests: only orgs with viol elec part
# all robustness tests are -- at first --
# only done with the main cross-sectional models and with the main specification
# if we were to fit all model specifications this would be
# at least: 4 methods * 3 inde. vars * 2 dep var * 3 hyptoheses = 72 models per robustness tests
library(specr)
library(fixest)
library(readr)
library(broom.mixed)
library(texreg)
library(dplyr)
library(fixest)

# Cross sectional linear FE models that use aggregate data (2000-2019)

source("funs_and_constants.R")

# load in data
epr_yearly_group_clusters_h1 <- read_csv("data/epro_data_h1_only_viol_elec.csv.gz")
epro_year_group_clusters <- read_csv("data/epro_data_h2_only_viol_elec.csv.gz")
epro_year_group_clusters_geo <- read_csv("data/epro_geo_data_h3_only_viol_elec.csv.gz")

mean_h1 <- epr_yearly_group_clusters_h1 %>% 
  filter(year >= 2000) %>%# only obs after 2000 because of EPR 2.0 coding
  group_by(gwgroupid, countries_gwid) %>%
  summarize(
    across(c(share_disagreement, disagreement, n_orgs, resource_and_agriculture, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic, nightlight_total_pc_log, n_tek_groups),
           ~ mean(.x, na.rm = T)),
    none_one_or_both_nats_agri_char = Mode(as.factor(none_one_or_both_nats_agri_char), na.rm = TRUE)
  )

mean_h2 <- epro_year_group_clusters %>% 
  filter(year >= 2000) %>%# only obs after 2000 because of EPR 2.0 coding
  group_by(gwgroupid, countries_gwid) %>%
  summarize(
    across(c(share_disagreement, disagreement, n_orgs, n_ed_religions, hhi_rel, religious_segments_bin, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic, nightlight_total_pc_log, n_tek_groups),
           ~ mean(.x, na.rm = T)))

mean_h3 <- epro_year_group_clusters_geo %>% 
  filter(year >= 2000) %>%# only obs after 2000 because of EPR 2.0 coding
  group_by(gwgroupid, countries_gwid) %>%
  summarize(
    across(c(share_disagreement, n_non_intersection_group_polygons, multiple_polygons_bin, spatial_hhi, disagreement, n_orgs, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic, nightlight_total_pc_log, n_tek_groups),
           ~ mean(.x, na.rm = T)))



# two main models per hypothesis

#######
## H1
#######


full_model_h1 <- feols(disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                       cluster = ~countries_gwid, data = mean_h1 %>% as.data.frame())

full_model_h1_n_resources <- feols(disagreement ~ none_one_or_both_nats_agri_char + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                                   cluster = ~countries_gwid, data = mean_h1)

screenreg(list(full_model_h1, full_model_h1_n_resources))


#######
### H2
#######

full_model_h2 <- feols(disagreement ~ religious_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                       cluster = ~countries_gwid, data = mean_h2)

# n rels as IV
full_model_h2_n_rels <- feols(disagreement ~ n_ed_religions + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                              cluster = ~countries_gwid, data = mean_h2)


full_model_h2_herf <- feols(disagreement ~ hhi_rel + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                            cluster = ~countries_gwid, data = mean_h2)

screenreg(list(full_model_h2, full_model_h2_n_rels, full_model_h2_herf))

#######
### H3
#######

full_model_h3 <- feols(disagreement ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                       cluster = ~countries_gwid, data = mean_h3)

full_model_h3_bin <- feols(disagreement ~ multiple_polygons_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                           cluster = ~countries_gwid, data = mean_h3)

full_model_h3_hhi <- feols(disagreement ~ spatial_hhi + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic + nightlight_total_pc_log + n_tek_groups | countries_gwid,
                           cluster = ~countries_gwid, data = mean_h3)

screenreg(list(full_model_h3, full_model_h3_bin, full_model_h3_hhi))


# make final table
texreg(list(full_model_h1, full_model_h1_n_resources, full_model_h2, full_model_h2_herf, full_model_h3, full_model_h3_hhi),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("H1" = 1:2, "H2" = 3:4, "H3" = 5:6),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/robustness_only_viol_elec.tex")

